Analizando datos speedy

(luego de la reunión del 5 de agosto de 2016)

Cambiando, decidimos dejar un poco de lado el geopotencial y ver más función corriente, u y T.

Variables: geop, streamf, T, U, V. Cosas para ver: * Estado medio * SD por latitud, mes * Mapas de las ondas * R^2 de las ondas * Cortes longitudinales de cada onda en distintas latitudes

Comparación con NCEP (pasar a la misma resolución usando interpolación bilineal)

Comparación con NCEP:

NCEP no brinda información de función corriente, para comparar uso geopotencial, U y T.

# Leo speedy y hago anomalías mensuales
m_speedy <- fread("SPEEDY/todo_m_speedy.csv")
m_speedy[, month := factor(month, levels = month.abb)]
##            lon     lat lev month         gh           u          v
##      1:   0.00 -87.159 925   Ene   736.8700 -0.05454388 -0.3241713
##      2:   3.75 -87.159 925   Ene   738.2792 -0.02213817 -0.3285298
##      3:   7.50 -87.159 925   Ene   739.8143  0.01149593 -0.3288868
##      4:  11.25 -87.159 925   Ene   741.4750  0.04583307 -0.3251453
##      5:  15.00 -87.159 925   Ene   743.2585  0.08036167 -0.3172835
##     ---                                                           
## 221180: 341.25  -1.856  30   Dic 23976.2846 -2.46443088  0.3448354
## 221181: 345.00  -1.856  30   Dic 23975.1531 -2.47587266  0.4439990
## 221182: 348.75  -1.856  30   Dic 23972.9128 -2.67345830  0.5403290
## 221183: 352.50  -1.856  30   Dic 23970.0830 -3.03758931  0.4697203
## 221184: 356.25  -1.856  30   Dic 23967.2616 -3.41034904  0.2759294
##                psi     temp
##      1: -0.1354184 267.7331
##      2: -0.1314101 267.6461
##      3: -0.1278538 267.5543
##      4: -0.1247234 267.4581
##      5: -0.1219879 267.3583
##     ---                    
## 221180: 15.9369396 206.9459
## 221181: 16.0817268 207.1822
## 221182: 16.2614006 207.4649
## 221183: 16.4383895 207.8276
## 221184: 16.5591864 208.0658
# setindex(speedy, lon, lat, lev, month, year)

vars <- c("gh", "u", "v", "psi", "temp")
avars <- paste0("a", vars)
allvars <- c(vars, avars)
sp_levs <- m_speedy[, unique(lev)]
# Anomalías mensuales 
m_speedy[, c(avars) := lapply(.SD, FUN = anom), 
         by = .(lat, lev, month),
         .SDcols = vars]
##            lon     lat lev month         gh           u          v
##      1:   0.00 -87.159 925   Ene   736.8700 -0.05454388 -0.3241713
##      2:   3.75 -87.159 925   Ene   738.2792 -0.02213817 -0.3285298
##      3:   7.50 -87.159 925   Ene   739.8143  0.01149593 -0.3288868
##      4:  11.25 -87.159 925   Ene   741.4750  0.04583307 -0.3251453
##      5:  15.00 -87.159 925   Ene   743.2585  0.08036167 -0.3172835
##     ---                                                           
## 221180: 341.25  -1.856  30   Dic 23976.2846 -2.46443088  0.3448354
## 221181: 345.00  -1.856  30   Dic 23975.1531 -2.47587266  0.4439990
## 221182: 348.75  -1.856  30   Dic 23972.9128 -2.67345830  0.5403290
## 221183: 352.50  -1.856  30   Dic 23970.0830 -3.03758931  0.4697203
## 221184: 356.25  -1.856  30   Dic 23967.2616 -3.41034904  0.2759294
##                psi     temp        agh         au          av
##      1: -0.1354184 267.7331  0.7799904 0.06318588 -0.36023388
##      2: -0.1314101 267.6461  2.1891375 0.09559159 -0.36459240
##      3: -0.1278538 267.5543  3.7242286 0.12922569 -0.36494939
##      4: -0.1247234 267.4581  5.3849729 0.16356284 -0.36120783
##      5: -0.1219879 267.3583  7.1684832 0.19809143 -0.35334602
##     ---                                                      
## 221180: 15.9369396 206.9459  5.0669393 1.14545083  0.14958942
## 221181: 16.0817268 207.1822  3.9354289 1.13400905  0.24875303
## 221182: 16.2614006 207.4649  1.6951294 0.93642341  0.34508308
## 221183: 16.4383895 207.8276 -1.1346883 0.57229240  0.27447440
## 221184: 16.5591864 208.0658 -3.9561076 0.19953267  0.08068342
##                  apsi       atemp
##      1: -0.0008850379  0.25983430
##      2:  0.0031231956  0.17283174
##      3:  0.0066795316  0.08099112
##      4:  0.0098099039 -0.01514434
##      5:  0.0125453975 -0.11492359
##     ---                          
## 221180:  2.6544129519 -0.83302260
## 221181:  2.7992001999 -0.59670922
## 221182:  2.9788739988 -0.31396132
## 221183:  3.1558628865  0.04871619
## 221184:  3.2766597577  0.28689470
setcolorder(m_speedy, c("lon", "lat", "lev", "month", allvars))
setorder(m_speedy, month, -lev, lat, lon)
m_ncep <- fread("NCEP/todo_m_ncep.csv")
m_ncep <- m_ncep[lat <= 0]
m_ncep[, month := factor(month, levels = month.abb)]
##            lon     lat  lev month           gh       temp           u
##      1:   0.00 -87.159 1000   Ene 3.194028e-01  -3.746012 -1.28944888
##      2:   3.75 -87.159 1000   Ene 2.009247e+00  -3.783916 -0.98486770
##      3:   7.50 -87.159 1000   Ene 3.678018e+00  -3.820472 -0.68162480
##      4:  11.25 -87.159 1000   Ene 5.273846e+00  -3.858975 -0.37856168
##      5:  15.00 -87.159 1000   Ene 6.822608e+00  -3.896709 -0.07922632
##     ---                                                              
## 470012: 341.25  -1.856   10   Dic 3.092088e+04 -42.738416 -7.69329709
## 470013: 345.00  -1.856   10   Dic 3.092209e+04 -42.620110 -7.63892995
## 470014: 348.75  -1.856   10   Dic 3.092438e+04 -42.696332 -7.64588785
## 470015: 352.50  -1.856   10   Dic 3.092421e+04 -42.711142 -7.58927508
## 470016: 356.25  -1.856   10   Dic 3.092315e+04 -42.988538 -7.17830892
##                   v
##      1: -3.37433997
##      2: -3.36687496
##      3: -3.33007242
##      4: -3.25966232
##      5: -3.16194629
##     ---            
## 470012: -0.17117253
## 470013: -0.02299725
## 470014: -0.08035137
## 470015: -0.24902995
## 470016: -0.65007687
m_ncep[, psi := NA]
##            lon     lat  lev month           gh       temp           u
##      1:   0.00 -87.159 1000   Ene 3.194028e-01  -3.746012 -1.28944888
##      2:   3.75 -87.159 1000   Ene 2.009247e+00  -3.783916 -0.98486770
##      3:   7.50 -87.159 1000   Ene 3.678018e+00  -3.820472 -0.68162480
##      4:  11.25 -87.159 1000   Ene 5.273846e+00  -3.858975 -0.37856168
##      5:  15.00 -87.159 1000   Ene 6.822608e+00  -3.896709 -0.07922632
##     ---                                                              
## 470012: 341.25  -1.856   10   Dic 3.092088e+04 -42.738416 -7.69329709
## 470013: 345.00  -1.856   10   Dic 3.092209e+04 -42.620110 -7.63892995
## 470014: 348.75  -1.856   10   Dic 3.092438e+04 -42.696332 -7.64588785
## 470015: 352.50  -1.856   10   Dic 3.092421e+04 -42.711142 -7.58927508
## 470016: 356.25  -1.856   10   Dic 3.092315e+04 -42.988538 -7.17830892
##                   v psi
##      1: -3.37433997  NA
##      2: -3.36687496  NA
##      3: -3.33007242  NA
##      4: -3.25966232  NA
##      5: -3.16194629  NA
##     ---                
## 470012: -0.17117253  NA
## 470013: -0.02299725  NA
## 470014: -0.08035137  NA
## 470015: -0.24902995  NA
## 470016: -0.65007687  NA
m_ncep[, temp := temp + 273.16]
##            lon     lat  lev month           gh     temp           u
##      1:   0.00 -87.159 1000   Ene 3.194028e-01 269.4140 -1.28944888
##      2:   3.75 -87.159 1000   Ene 2.009247e+00 269.3761 -0.98486770
##      3:   7.50 -87.159 1000   Ene 3.678018e+00 269.3395 -0.68162480
##      4:  11.25 -87.159 1000   Ene 5.273846e+00 269.3010 -0.37856168
##      5:  15.00 -87.159 1000   Ene 6.822608e+00 269.2633 -0.07922632
##     ---                                                            
## 470012: 341.25  -1.856   10   Dic 3.092088e+04 230.4216 -7.69329709
## 470013: 345.00  -1.856   10   Dic 3.092209e+04 230.5399 -7.63892995
## 470014: 348.75  -1.856   10   Dic 3.092438e+04 230.4637 -7.64588785
## 470015: 352.50  -1.856   10   Dic 3.092421e+04 230.4489 -7.58927508
## 470016: 356.25  -1.856   10   Dic 3.092315e+04 230.1715 -7.17830892
##                   v psi
##      1: -3.37433997  NA
##      2: -3.36687496  NA
##      3: -3.33007242  NA
##      4: -3.25966232  NA
##      5: -3.16194629  NA
##     ---                
## 470012: -0.17117253  NA
## 470013: -0.02299725  NA
## 470014: -0.08035137  NA
## 470015: -0.24902995  NA
## 470016: -0.65007687  NA
m_ncep[,  c(avars) := lapply(.SD, FUN = anom), by = .(lat, lev, month), 
            .SDcols = vars]
##            lon     lat  lev month           gh     temp           u
##      1:   0.00 -87.159 1000   Ene 3.194028e-01 269.4140 -1.28944888
##      2:   3.75 -87.159 1000   Ene 2.009247e+00 269.3761 -0.98486770
##      3:   7.50 -87.159 1000   Ene 3.678018e+00 269.3395 -0.68162480
##      4:  11.25 -87.159 1000   Ene 5.273846e+00 269.3010 -0.37856168
##      5:  15.00 -87.159 1000   Ene 6.822608e+00 269.2633 -0.07922632
##     ---                                                            
## 470012: 341.25  -1.856   10   Dic 3.092088e+04 230.4216 -7.69329709
## 470013: 345.00  -1.856   10   Dic 3.092209e+04 230.5399 -7.63892995
## 470014: 348.75  -1.856   10   Dic 3.092438e+04 230.4637 -7.64588785
## 470015: 352.50  -1.856   10   Dic 3.092421e+04 230.4489 -7.58927508
## 470016: 356.25  -1.856   10   Dic 3.092315e+04 230.1715 -7.17830892
##                   v psi       agh         au          av apsi       atemp
##      1: -3.37433997  NA 10.054246 -0.0262046 -4.17046447   NA -0.05855587
##      2: -3.36687496  NA 11.744090  0.2783766 -4.16299946   NA -0.09645984
##      3: -3.33007242  NA 13.412862  0.5816195 -4.12619692   NA -0.13301612
##      4: -3.25966232  NA 15.008689  0.8846826 -4.05578682   NA -0.17151928
##      5: -3.16194629  NA 16.557452  1.1840180 -3.95807079   NA -0.20925289
##     ---                                                                  
## 470012: -0.17117253  NA  3.675745  2.1369449 -0.04305793   NA  0.19750639
## 470013: -0.02299725  NA  4.888170  2.1913120  0.10511735   NA  0.31581262
## 470014: -0.08035137  NA  7.173449  2.1843541  0.04776323   NA  0.23959038
## 470015: -0.24902995  NA  7.003035  2.2409669 -0.12091535   NA  0.22478026
## 470016: -0.65007687  NA  5.947643  2.6519330 -0.52196226   NA -0.05261632
np_levs <- m_ncep[, unique(lev)]
setcolorder(m_ncep, c("lon", "lat", "lev", "month", allvars))
setorder(m_ncep, month, -lev, lat, lon)
m_speedy.long <- melt(m_speedy, id.vars = c("lon", "lat", "lev", "month"))
m_speedy.interp <- m_speedy.long[, 
                    approx(x = lev, y = value,
                              xout = np_levs),
                    by = .(lon, lat, month, variable)]
colnames(m_speedy.interp)[5:6] <- c("lev", "sp")
m_speedy <- m_speedy.interp

Correlación lineal para cada nivel, mes y variable, separando en parte total y parte asimétrica.

m_ncep.long <- melt(m_ncep, id.vars = c("lon", "lat", "lev", "month"), value.name = "nc")
sp_nc <- merge(m_speedy.interp, m_ncep.long)

remove(m_speedy.interp)

cors <- sp_nc[lev %in% sp_levs, .(Correlacion = cor(sp, nc)), by = .(lev, month, variable)]
levs <- unique(cors$lev)

corsa <- cors[grepl("a", variable)]$Correlacion
cors <- cors[!grepl("a", variable)]
cors[, Asimetrica := corsa]
##      lev month variable Correlacion Asimetrica
##   1:  30   Ene       gh  -0.8679965  0.6589105
##   2: 100   Ene       gh   0.9841290  0.7028101
##   3: 200   Ene       gh   0.9856019  0.7279014
##   4: 300   Ene       gh   0.9838896  0.6130320
##   5: 500   Ene       gh   0.9801555  0.4849896
##  ---                                          
## 476: 300   Dic     temp   0.9751213  0.5669539
## 477: 500   Dic     temp   0.9929736  0.5298177
## 478: 700   Dic     temp   0.9819914  0.3370901
## 479: 850   Dic     temp   0.9843473  0.6985642
## 480: 925   Dic     temp   0.9895365  0.8309316
colnames(cors) <- c("lev", "month", "variable", "Total", "Asimetrica")
cors <- melt(cors, id.vars = c("lev", "month", "variable"), variable.name = "Parte", 
              value.name = "Correlacion")
ggplot(cors, aes(lev, Correlacion, color = Parte)) +
  geom_hline(yintercept = 0, color = "gray45") +
  geom_line() + facet_grid(month~variable) +
  scale_color_brewer(type = "qual", name = "Variable") +
  coord_flip() + scale_x_continuous(trans = "reverselog", breaks = levs) +
  scale_y_continuous(limits = c(-1, 1), minor_breaks = NULL) +
  xlab("Nivel") +  theme_bw() +
  theme(legend.position = "bottom", legend.key.height = unit(5, "points")) +
  ggtitle("Correlación lineal entre Speedy y NCEP ")

## Note: no visible binding for global variable 'xend' 
## Note: no visible binding for global variable 'yend' 
## Note: no visible binding for global variable 'x' 
## Note: no visible binding for global variable 'y' 
## Note: no visible binding for global variable 'R' 
## Note: no visible binding for global variable 'R' 
## Note: no visible binding for global variable 'R' 
## Note: no visible binding for global variable 'R' 
## Note: no visible binding for global variable 'R' 
## Note: no visible binding for global variable 'R' 
## Note: no visible binding for global variable 'R' 
## Note: no visible binding for global variable 'R' 
## Note: no visible binding for global variable 'R' 
## Note: no visible binding for global variable 'R' 
## Note: no visible binding for global variable 'R' 
## Note: no visible binding for global variable 'R' 
## Note: no visible binding for global variable 'R' 
## Note: no visible binding for global variable 'R' 
## Note: no visible binding for global variable 'R' 
## Note: no visible binding for global variable 'R' 
## Note: no visible binding for global variable 'key.row' 
## Note: no visible binding for global variable 'key.col' 
## Note: no visible binding for global variable 'label.row' 
## Note: no visible binding for global variable 'label.col' 
## Note: no visible binding for global variable 'key.row' 
## Note: no visible binding for global variable 'key.col' 
## Note: no visible binding for global variable 'label.row' 
## Note: no visible binding for global variable 'label.col' 
## Note: no visible binding for global variable 'key.row' 
## Note: no visible binding for global variable 'key.col' 
## Note: no visible binding for global variable 'label.row' 
## Note: no visible binding for global variable 'label.col' 
## Note: no visible binding for global variable 'key.row' 
## Note: no visible binding for global variable 'key.col' 
## Note: no visible binding for global variable 'label.row' 
## Note: no visible binding for global variable 'label.col'
interp.levs <- exp(seq(log(925), log(30), length.out = 40))

cors.interp <- cors[!is.na(Correlacion), 
                    interp.dt(as.numeric(month), lev, Correlacion, linear = T,
                              yo = interp.levs),
                    by = .(Parte, variable)]
## Note: no visible binding for global variable 'x' 
## Note: no visible binding for global variable 'y' 
## Note: no visible binding for global variable 'z'
colnames(cors.interp) <- c("Parte", "variable", "month", "lev", "Correlacion")

cor_plot <- ggplot(cors.interp, aes(month, lev)) +
  geom_raster(aes(fill = Correlacion)) +
  geom_contour(aes(z = Correlacion, color = ..level..), 
               binwidth = 0.2) +
  facet_wrap(variable~Parte, scales = "free", ncol = 2) + 
  scale_y_continuous(trans = "reverselog", breaks = levs) +
  scale_x_continuous(labels = month.abb, breaks = 1:12) +
  scale_fill_distiller(type = "div", palette = "RdBu", direction = 1) +
  scale_color_gradient(low = "black", high = "black") +
  theme_minimal()
  # scale_fill_gradient2(low = "#67001f", high = "#053061")
cor_plot <- direct.label(cor_plot, method = "bottom.pieces")
  • Altura geopotencial (gh) Para el caso del campo total, la correlación del campo es buena (>0.8)en casi todos los niveles y meses, excepto en 30hPa durante verano donde los campos ¡está anticorrelacionados! La parte asimétrica zonal muestra valores menores, indicando que gran parte de la correlación del campo total se debe a la capacidad del modelo de reproducir el gradiente meridional. Sin embargo, se siguen obteniendo correlaciones >0.6 en casi todos los niveles y meses. Se observa un mínimo relativo en 500hPa donde en se tienen correlaciones menores durante casi todo el año y uno en niveles altos centrado en septiembre, donde la correlación llega a ser nula.

  • Viento zonal (U) Las correlaciones con el campo total son >=0.8 en todo el año y todos los niveles, sin embargo, la parte asimétrica muestra correlaciones mucho más baja con un máximo de ~0.6 en 925hPa. Esto indica que el modelo resuelve correctamente la estructura media del Jet, pero no sus variaciones zonales.

  • Viento meridional (V) Los campos de correlación son prácticamente idénticos entre parte total y parte asimétrica. Ésta muestra un patrón de bajas correlaciones en general, con anticorrelaciones en niveles altos (entre 300 y 200 hPa) durante todo el año

  • Temperatura (T) La correlación con el campo total muestra una estructura similar que la altura geopotencial, con una excelente correlación en todos los meses para niveles mayores a 200hPa, pero anticorrelacionado en niveles altos en todos los meses salvo en invierno. La parte asimétrica muestra correlaciones bajas en todos los niveles salvo en 925hPa.

Parte simétrica

Altura geopotencial

sp_nc_mean <- sp_nc[, lapply(.SD, mean), by = .(lat, lev, month, variable), 
                    .SDcols = c("sp", "nc")]
sp_nc_mean.long <- melt(sp_nc_mean, id.vars = c("lat", "lev", "month", "variable"), variable.name = "modelo")


ggplot(sp_nc_mean.long[variable == "gh"], aes(lat, value/1000)) + 
  geom_line(aes(color = modelo, group = interaction(modelo, as.factor(lev)), 
                linetype = as.factor(lev))) +
  scale_linetype_manual(values = rep(c(1,2), length.out = 17)) +
  guides(linetype = FALSE) +
  geom_text(data = sp_nc_mean.long[month == "Dic" & lat %~% -50 & modelo == "sp" & variable == "gh"], aes(label = lev), nudge_y = 0) +
  facet_grid(~month) + 
  theme_elio + 
  scale_color_brewer(type = "qual", palette = 6, labels = c("Speedy", "NCEP"),
                     name = "Modelo") +
  scale_x_reverse() +
  ylab("Media zonal de altura geopotencial (kmgp)") + xlab("Latitud")

interp.levs <- exp(seq(log(1000), log(10), length.out = 60))
sp_nc_corte <- sp_nc_mean.long[!is.na(value), 
                               interp.dt(lat, lev, value, linear = T,
                                         yo = interp.levs, nx = 60),
                               by = .(month, variable, modelo)]
colnames(sp_nc_corte)[4:6] <- c("lat", "lev", "value")
sp_nc_corte.wide <- data.table::dcast(sp_nc_corte, month + variable + lat + lev ~ modelo)

corte_latlev_zonal <- function(var, lineas_width, title, scale = "seq") {
  M <- sp_nc_corte[variable == var, max(value, na.rm = T)]
  m <- sp_nc_corte[variable == var, min(value,  na.rm = T)]
  m <- floor(m/lineas_width)*lineas_width
  M <- ceiling(M/lineas_width)*lineas_width
  
  g <- ggplot(sp_nc_corte.wide[variable == var], aes(lat, lev))
   
  # lineas_width <- 5000
  if (scale == "div") {
    side <- max(-m, M)
    lineas = seq(-side, side, by = lineas_width)
    div_pal = colorRampPalette(brewer.pal(name = "RdBu", n = 11))
    colors <- div_pal(length(lineas) - 1)
    colors <- colors[length(colors):1]
    inicio <- 0
  g <- g + geom_raster(aes(fill = nc)) + 
    scale_fill_manual(values = colors, drop = FALSE, name = "NCEP") 
  } else {
    lineas = seq(m, M, by = lineas_width)
    inicio <- lineas[1]
    g <- g + geom_raster(aes(fill = cut_width(nc, lineas_width, boundary = inicio))) +
      scale_fill_viridis(name = "NCEP", option = "plasma",
                       direction = 1, discrete = T, drop = FALSE,
                       label = c(lineas, lineas[length(lineas)] + lineas_width))
  }
  
 g <- g +  
    geom_contour(aes(z = sp), alpha = 0.8, color = "black", breaks = lineas) +
    geom_dl(aes(z = sp, label = ..level..), method = "top.pieces", 
            stat = "contour", color = "black", breaks = jumpby(lineas, 2)) +
    facet_wrap(~month, ncol = 3) +
    scale_y_continuous(trans = "reverselog", breaks = levs) +
    scale_x_reverse() +
    theme_minimal() + 
    ggtitle(paste0(title, "\n lineas = ", lineas_width)) +
    geom_hline(yintercept = 30, color = "gray45", linetype = 2)
}

g <- corte_latlev_zonal(var = "gh", lineas_width = 5000, title = "Altura geopotencial media (mgp)")
## Note: no visible binding for global variable 'variable' 
## Note: no visible binding for global variable 'value' 
## Note: no visible binding for global variable 'variable' 
## Note: no visible binding for global variable 'value' 
## Note: no visible binding for global variable 'variable' 
## Note: no visible binding for global variable 'lat' 
## Note: no visible binding for global variable 'lev' 
## Note: no visible binding for global variable 'nc' 
## Note: no visible binding for global variable 'nc' 
## Note: no visible binding for global variable 'sp' 
## Note: no visible binding for global variable 'sp' 
## Note: no visible binding for global variable '..level..'
g

Comparando la parte zonalmente simétrica de la altura geopotencial se hace evidente la razón de la anticorrelación en niveles altos durante verano. En niveles bajos y medios, el gradiente meridional de geopotencial es positivo (más alturas hacia el ecuador) durane todo el año, con una amplitid que aumenta con la altura. En niveles altos, sin embargo, durante verano NCEP muestra que éste se invierte, mostrando mayores alturas en el altas latitudes. Speedy, por su parte, no logra capturar este comportamiento y sigue repitiendo la estructura básica de nieles inferiores. Esto puede deberse a que Speedy tiene una estratósfera en un nivel más alto que NCEP en los meses de verano.

Más allá de esta salvedad, se ve que todos los niveles existe una buena concordancia entre ambos modelos. La excepción es 50 y 70 hPa, donde Speedy sobreestima la altura geopotencial en todos los meses y latitudes considerablemente.

Temperatura

En la temperatura, se observa que la anticorrelación que se observa en 200 hPa tiene una explicación similar a la de la altura geopotencial. En niveles bajos, el gradiente de temperatura zonal es positivo (mayores temperaturas en el ecuador) mientras que en niveles altos, por encima de 100hPa, este gradiente de invierte. La región intermedia, donde el gradiente es nulo se encuentra cerca de los 200hPa. En los meses de verano para Speedy ésta está lieramente más alta que para NCEP, y ese defasaje es lo que da lugar a la anticorrelación.

Esta diferencia puede indicar que Speedy no captura algún proceso estratosférico. ¿Será que lo que se observa es producto del calentamiento por la producción de Ozono? Hay que investigar cómo Speedy trata el proceso de ozonogénesis, ya que podría ser que no tiene en cuenta el calor extra o lo subestima.

Analizando la parte simétrica de la temperatura, es evidente que Speedy subestima la temperatura significativamente, con diferencias de casi 60K en bajas latitudes en 100hPa. Por otro lado, también sbreestima el ciclo anual en altas latitudes; tanto que en verano invierte el gradiente de temperatura, dando origen a la anticorrelación observada. (En otras palabras, es un desastre)

Observando el perfil en altura.. ¿estratósfera? En -90 speedy parece que ve la estratósfera en en 30hPa mientras que NCEP la pone mucho más abajo, en 300hPa en invierno. En trópico y ecuador no se aprecia la inversión estratosférica.

Speedy sobreestima el gradiente de temperatura. ¿Esto significa más inestabilidad?

plotlevs <- c(850, 500, 300, 200, 100, 30)
ggplot(sp_nc_mean.long[variable == "temp" & lev %in% plotlevs], aes(lat, value)) + 
  geom_line(aes(color = modelo, group = interaction(modelo, as.factor(lev)))) +
  facet_grid(lev~month) + 
  theme_minimal() + 
  scale_color_brewer(type = "qual", palette = 6, labels = c("Speedy", "NCEP"),
                     name = "Modelo") +
  scale_y_continuous() +
  scale_x_reverse() +
  ylab("Media zonal de la temperatura (K)") + xlab("Latitud")

g <- corte_latlev_zonal(var = "temp", lineas_width = 10, title = "Temperatura media (K)")
g

lats <- sp_nc_mean.long[, unique(lat)][c(1, 6, 12, 18, 24)]

ggplot(sp_nc_mean.long[variable == "temp" & lat %in% lats], aes(lev, value - 273.16)) +
  geom_line(aes(color = modelo, group = interaction(modelo, as.factor(lat)))) +
  facet_grid(lat~month) +
  coord_flip() +
  theme_minimal() +
  scale_color_brewer(type = "qual", palette = 6, labels = c("Speedy", "NCEP"),
                     name = "Modelo") +
  scale_x_continuous(trans = "reverselog", breaks = levs) +
  scale_y_continuous(breaks = seq(-90, 30, by = 40)) +
  xlab("Nivel (hPa)") + ylab("Temperatura  (C)")

Viento zonal

# g <- corte_latlev_zonal(var = "u", lineas_width = 10, title = "Viento zonal (m/s)", scale = "div")
# 
# div_pal = colorRampPalette(brewer.pal(name = "RdBu", n = 11))
# g + discrete_scale(values = div_pal(20), name = "NCEP", drop = FALSE)
# 
#   
plot_var = "u"
M <- ceiling(sp_nc_corte[variable == plot_var, max(value, na.rm = T)])
m <- floor(sp_nc_corte[variable == plot_var, min(value,  na.rm = T)])
lineas_width <- 10
m <- floor(m/lineas_width)*lineas_width
lineas = seq(m, M, by = lineas_width)

ggplot(sp_nc_corte.wide[variable == plot_var], aes(lat, lev)) +
  geom_raster(aes(fill = nc)) +
  geom_contour(aes(z = sp, linetype = as.factor(-sign(..level..))),
               alpha = 0.8, color = "black", breaks = lineas) +
  geom_dl(aes(z = sp, label = ..level..), method = "top.pieces", 
          stat = "contour", color = "black", breaks = jumpby(lineas, 2)) +
  
  geom_contour(aes(z = nc), 
               alpha = 0.8, color = "grey33", breaks = lineas) +
  geom_dl(aes(z = nc, label = ..level..), method = "bottom.pieces", 
          stat = "contour", color = "grey33", breaks = jumpby(lineas, 2)) +
  scale_linetype_manual(values = c(1, 2, 3)) +
  facet_wrap(~month, ncol = 3) +
  scale_y_continuous(trans = "reverselog", breaks = levs) +
  scale_x_reverse() +
  scale_fill_gradient2(name = "NCEP", high = muted("red"), low = muted("blue")) +
  theme_minimal() + ggtitle(paste0("Viento zonal (m/s) \n lineas = ", lineas_width))

Por su parte, el viento zonal ambos modelos muestran una estructura similar con el jet subtropical bien definido aunque Speedy lo muestra más intenso y ligeramente corrido hacia el ecuador y en niveles más altos que NCEP durante los meses de verano. Durante el invierno, NCEP muestra el jet subpolar, pero Speedy no logra representarlo por falta de niveles verticales, aunque la estructura por debajo del jet se muestra similar aunque subestimando la magnitud.

Speedy tampoco muestra los vientos del este en niveles altos que dominan bajas latitudes en invierno y casi todas las latitudes en verano. Tiene un mínimo y valores negativos, pero muy subestimados.

ggplot(sp_nc_mean.long[variable == "u" & lev %in% plotlevs], aes(lat, value)) + 
  geom_line(aes(color = modelo, group = interaction(modelo, as.factor(lev)))) +
  facet_grid(lev~month) + 
  theme_minimal() + 
  scale_color_brewer(type = "qual", palette = 6, labels = c("Speedy", "NCEP"),
                     name = "Modelo") +
  scale_y_continuous(breaks = seq(-10, 60, by = 20)) +
  scale_x_reverse() +
  ylab("Media zonal del viento zonal (m/s)") + xlab("Latitud") +
  theme_elio

Viento meridional

plot_var = "v"
M <- ceiling(sp_nc_corte[variable == plot_var, max(value, na.rm = T)])
m <- floor(sp_nc_corte[variable == plot_var, min(value,  na.rm = T)])
lineas_width <- .5
m <- floor(m/lineas_width)*lineas_width
lineas = seq(lineas_width, M, by = lineas_width)

ggplot(sp_nc_corte.wide[variable == plot_var], aes(lat, lev)) +
  geom_raster(aes(fill = nc)) +
  geom_contour(aes(z = sp, linetype = as.factor(sign(..level..))),
               alpha = 0.9, color = "black", breaks = fill_neg(lineas)) +
  geom_dl(aes(z = sp, label = ..level..), method = "top.pieces", 
          stat = "contour", color = "black", breaks = jumpby(fill_neg(lineas), 2)) +
  
  geom_contour(aes(z = nc), 
               alpha = 0.6, color = "grey33", breaks = fill_neg(lineas)) +
  geom_dl(aes(z = nc, label = ..level..), method = "bottom.pieces", 
          stat = "contour", color = "grey33", breaks = jumpby(fill_neg(lineas), 2)) +
  scale_linetype_manual(values = c(2, 1, 1), guide = "none") +
  facet_wrap(~month, ncol = 3) +
  scale_y_continuous(trans = "reverselog", breaks = levs) +
  scale_x_reverse() +
  scale_fill_gradient2(name = "NCEP", high = muted("red"), low = muted("blue")) +
  theme_elio + ggtitle(paste0("Viento zonal (m/s) \n lineas = ", lineas_width))

plot_var = "v"
M <- ceiling(sp_nc_corte[variable == plot_var, max(value, na.rm = T)])
m <- floor(sp_nc_corte[variable == plot_var, min(value,  na.rm = T)])
lineas_width <- .5
m <- floor(m/lineas_width)*lineas_width
lineas = seq(lineas_width, M, by = lineas_width)

ggplot(sp_nc_corte.wide[variable == plot_var & lev >= 500], aes(lat, lev)) +
  geom_raster(aes(fill = nc)) +
  geom_contour(aes(z = sp, linetype = as.factor(sign(..level..))),
               alpha = 0.9, color = "black", breaks = fill_neg(lineas)) +
  geom_dl(aes(z = sp, label = ..level..), method = "top.pieces", 
          stat = "contour", color = "black", breaks = jumpby(fill_neg(lineas), 2)) +
  
  geom_contour(aes(z = nc), 
               alpha = 0.6, color = "grey33", breaks = fill_neg(lineas)) +
  geom_dl(aes(z = nc, label = ..level..), method = "bottom.pieces", 
          stat = "contour", color = "grey33", breaks = jumpby(fill_neg(lineas), 2)) +
  scale_linetype_manual(values = c(2, 1, 1), guide = "none") +
  facet_wrap(~month, ncol = 3) +
  scale_y_continuous(trans = "reverselog", breaks = levs) +
  scale_x_reverse() +
  scale_fill_gradient2(name = "NCEP", high = muted("red"), low = muted("blue")) +
  theme_elio + ggtitle(paste0("Viento zonal (m/s) \n lineas = ", lineas_width))

En el viento meridional, Speedy parece capturar la estructura general de la circulación pero subestima las magnitudes. En bajas latitudes, hay un dipolo entre niveles bajos y altos que alterna entre invierno y verano; se trata de la parte convergente en superficie y divergente en altura de la ITCZ, que se mueve hacia el hemisferio de verano.

En altas latitudes, en superficie hay máximos de viento del sur debido a los vientos catabáticos de la antártida. Nuevamente, Speedy captura esta característica, pero con una considerable subestimación.

ggplot(sp_nc_mean.long[variable == "v" & lev %in% plotlevs], aes(lat, value)) + 
  geom_line(aes(color = modelo, group = interaction(modelo, as.factor(lev)))) +
  facet_grid(lev~month, scales = "free_y") + 
  theme_elio + 
  scale_color_brewer(type = "qual", palette = 6, labels = c("Speedy", "NCEP"),
                     name = "Modelo") +
  scale_y_continuous() +
  scale_x_reverse() +
  ylab("Media zonal de la temperatura (K)") + xlab("Latitud")

Parte Antisimétrica

Altura geopotencial

Campos por estación y nivel:

sp_nc[, estacion := asign_season(month)]
##             lon     lat month variable  lev            sp            nc
##       1:   0.00 -87.159   Ene       gh   10            NA 31955.3946220
##       2:   0.00 -87.159   Ene       gh   20            NA 27002.9561155
##       3:   0.00 -87.159   Ene       gh   30 23577.7435547 24162.6416817
##       4:   0.00 -87.159   Ene       gh   50 21357.7593099 20622.8639138
##       5:   0.00 -87.159   Ene       gh   70 19137.7750651 18308.5150735
##      ---                                                               
## 4700156: 356.25  -1.856   Dic    atemp  600    -0.4605941    -0.5462326
## 4700157: 356.25  -1.856   Dic    atemp  700    -1.5416182    -0.5649992
## 4700158: 356.25  -1.856   Dic    atemp  850    -1.3644017    -0.4677048
## 4700159: 356.25  -1.856   Dic    atemp  925    -1.8917754    -0.4087908
## 4700160: 356.25  -1.856   Dic    atemp 1000            NA    -0.8071430
##          estacion
##       1:   Verano
##       2:   Verano
##       3:   Verano
##       4:   Verano
##       5:   Verano
##      ---         
## 4700156:   Verano
## 4700157:   Verano
## 4700158:   Verano
## 4700159:   Verano
## 4700160:   Verano
g_data <- rep_lon(sp_nc[variable == "agh", .(sp = mean(sp), nc = mean(nc)),
                        by = .(lat, lon, lev, estacion)])
## Note: no visible binding for global variable 'lon' 
## Note: no visible binding for global variable 'lon'
maplevs <- plotlevs[plotlevs != 850]
m_data <- g_data[lev %in% maplevs]
M <- m_data[, (max(c(abs(min(nc)), max(nc))))]
M <- ceiling(M)
# cuts <- seq(-M, M, length.out = 10)
# lineas <- floor(abs(cuts[1] - cuts[2]))
# cuts <- cuts[cuts != 0]
s <- 0.5
linea <- 60
lineas <- seq(linea, M, by = linea)
# m_data <- m_data[lev == 30 & estacion == "Primavera"]

ggplot(m_data, aes(lon, lat)) +
  # geom_tile(aes(fill = nc)) +
  facet_grid(lev~estacion) +
  HS3 +
  geom_contour(aes(z = sp, linetype = as.factor(sign(..level..))),
               breaks = fill_neg(lineas), color = "red4", size = s, alpha = 0.8) +
  # geom_dl(aes(z = sp, label = ..level..), method = "top.pieces", stat = "contour", breaks = fill_neg(jumpby(lineas, 2)), color = "red4") + 
  
  geom_contour(aes(z = nc, linetype = as.factor(sign(..level..))),
               breaks = fill_neg(lineas), color = "blue4", size = s,  alpha = 0.8) +
  # geom_dl(aes(z = nc, label = ..level..), method = "bottom.pieces", stat = "contour", breaks = fill_neg(jumpby(lineas, 2)), color = "blue4") + 
  
  scale_fill_gradient2(high = muted("red"), low = muted("blue"), 
                       name = "NCEP") +
  scale_linetype_manual(values = c(2, 1, 1)) +
  coord_map("stereographic", orientation = c(-90,0, 60),
            ylim = c(-90, -20)) +
  scale_x_continuous(breaks = seq(0, 330, by = 30)) +
  # annotate(geom = "text", x = seq(0, 330, by = 30), y = -5, 
  #          label = seq(0, 330, by = 30), color = "gray45") +
  # scale_y_continuous(breaks = seq(0+15, -90+15, by = -15)) +
  # annotate(geom = "text", x = 0, y = seq(0+15, -90+15, by = -15), 
  #          label = seq(0+15, -90+15, by = -15), color = "gray45") +
  theme_elio + 
  theme(legend.position = "bottom", legend.key.height = unit(5, "points"),
        axis.title = element_blank(), axis.text = element_blank()) +
  ggtitle(paste0("Z* para NCEP (azul) y SPEEDY (rojo) \n lineas = ", linea))

Corte en -65º

ggplot(g_data[lat %~% -65], aes(lon, lev)) +
  geom_contour(aes(z = sp, linetype = as.factor(sign(..level..))),
               breaks = fill_neg(lineas), color = "red4", size = s,  alpha = 0.8) +
  geom_contour(aes(z = nc, linetype = as.factor(sign(..level..))), 
               breaks = fill_neg(lineas), color = "blue4", size = s,  alpha = 0.8) +
  scale_linetype_manual(values = c(2, 1, 1), guide = "none") +
  facet_grid(~estacion) +
  scale_y_continuous(trans = "reverselog", breaks = levs) +
  theme_bw() +
  ggtitle(paste0("Z* para NCEP (azul) y SPEEDY (rojo) en -65º \n lineas = ", linea)) 

Un pequeño experimento no del todo exitoso para visualizar el corte:

ggplot(g_data[lat %~% -65], aes(lon, lev)) +
  geom_contour(aes(z = sp, linetype = as.factor(sign(..level..))),
               breaks = fill_neg(lineas), color = "red4", size = s,  alpha = 0.8) +
  geom_contour(aes(z = nc, linetype = as.factor(sign(..level..))), 
               breaks = fill_neg(lineas), color = "blue4", size = s,  alpha = 0.8) +
  scale_linetype_manual(values = c(2, 1, 1), guide = "none") +
  facet_grid(~estacion) +
  scale_y_continuous(trans = "reverselog", breaks = levs, limits = c(15000, 10)) +
  scale_x_continuous(breaks = seq(60, 360, by = 60)) +
  theme_bw() +
  ggtitle(paste0("Z* para NCEP (azul) y SPEEDY (rojo) en -65º \n lineas = ", linea)) +
   coord_polar(start = -300*pi/180, direction = 1) + # annotate(geom = "line", x = c(0, 360), y = rep(925, 2)) +
  annotate(geom = "rect", xmax = 360, xmin = 0, ymax = 15000, ymin = 925, fill = "white")

## Note: no visible global function definition for 'scale_transform'

Desvío estándar de anomalía de geopotencial.

sp_nc_sd <- sp_nc[, .(sp = sd(sp), nc = sd(nc)), 
                                     by = .(lat, month, lev, variable)]
g_data <- sp_nc_sd[lev %in% plotlevs & variable == "agh"]

M <- g_data[, (max(c(abs(min(nc)), max(nc))))]
M <- ceiling(M)
s <- 0.7
linea <- 40
lineas <- seq(linea, M, by = linea)

ggplot(g_data, aes(as.numeric(month), lat)) +
  facet_grid(lev~.) +
  geom_contour(aes(z = sp), breaks = fill_neg(lineas), 
               color = "red4", size = s,  alpha = 0.8) +
  geom_dl(aes(label = ..level.., z = sp), color = "red4", 
          stat = "contour", method = "bottom.pieces", 
          breaks = jumpby(fill_neg(lineas), 2)) +
  
  geom_contour(aes(z = nc), breaks = fill_neg(lineas), 
               color = "blue4", size = s,  alpha = 0.8) +
  geom_dl(aes(label = ..level.., z = nc), color = "blue4", 
          stat = "contour", method = "top.pieces", 
          breaks = jumpby(fill_neg(lineas), 2)) + 
  theme_elio +
  scale_x_continuous(breaks = 1:12, labels = month.abb)

Viento zonal.

g_data <- rep_lon(sp_nc[variable == "au", .(sp = mean(sp), nc = mean(nc)),
                        by = .(lat, lon, lev, estacion)])
maplevs <- plotlevs[plotlevs != 850]
m_data <- g_data[lev %in% maplevs]
M <- m_data[, (max(c(abs(min(nc)), max(nc))))]
M <- ceiling(M)

s <- 0.3
linea <- 5
lineas <- seq(linea, M, by = linea)

ggplot(m_data, aes(lon, lat)) +
  # geom_tile(aes(fill = nc)) +
  facet_grid(lev~estacion) +
  HS3+
  geom_contour(aes(z = sp, linetype = as.factor(sign(..level..))),
               breaks = fill_neg(lineas), color = "red4", size = s, alpha = 0.8) +
  # geom_dl(aes(z = sp, label = ..level..), method = "top.pieces", stat = "contour", breaks = fill_neg(jumpby(lineas, 2)), color = "red4") + 
  
  geom_contour(aes(z = nc, linetype = as.factor(sign(..level..))),
               breaks = fill_neg(lineas), color = "blue4", size = s,  alpha = 0.8) +
  # geom_dl(aes(z = nc, label = ..level..), method = "bottom.pieces", stat = "contour", breaks = fill_neg(jumpby(lineas, 2)), color = "blue4") + 
  
  scale_fill_gradient2(high = muted("red"), low = muted("blue"), 
                       name = "NCEP") +
  scale_linetype_manual(values = c(2, 1, 1), guide = "none") +
  coord_map("stereographic", orientation = c(-90,0, 60),
            ylim = c(-90, 0)) +
  scale_x_continuous(breaks = seq(0, 330, by = 30)) +
  # annotate(geom = "text", x = seq(0, 330, by = 30), y = -5, 
  #          label = seq(0, 330, by = 30), color = "gray45") +
  # scale_y_continuous(breaks = seq(0+15, -90+15, by = -15)) +
  # annotate(geom = "text", x = 0, y = seq(0+15, -90+15, by = -15), 
  #          label = seq(0+15, -90+15, by = -15), color = "gray45") +
  theme_bw() + 
  theme(legend.position = "bottom", legend.key.height = unit(5, "points"),
        axis.title = element_blank(), axis.text = element_blank()) +
  ggtitle(paste0("U* para NCEP (azul) y SPEEDY (rojo) \n lineas = ", linea))

linea <- 2
lineas <- seq(linea, M, by = linea)

ggplot(g_data[lat %~% -65], aes(lon, lev)) +
  geom_contour(aes(z = sp, linetype = as.factor(sign(..level..))),
               breaks = fill_neg(lineas), color = "red4", size = s) +
  geom_contour(aes(z = nc, linetype = as.factor(sign(..level..))), 
               breaks = fill_neg(lineas), color = "blue4", size = s) +
  scale_linetype_manual(values = c(2, 1, 1), guide = "none") +
  facet_grid(~estacion) +
  scale_y_continuous(trans = "reverselog", breaks = levs) +
  theme_bw() +
  ggtitle(paste0("U* para NCEP (azul) y SPEEDY (rojo) en -65º \n lineas = ", linea))

g_data <- sp_nc_sd[lev %in% plotlevs & variable == "au"]

M <- g_data[, (max(c(abs(min(nc)), max(nc))))]
M <- ceiling(M)
s <- 0.7
linea <- 4.5
lineas <- seq(linea, M, by = linea)

ggplot(g_data, aes(as.numeric(month), lat)) +
  facet_grid(lev~.) +
  geom_contour(aes(z = sp), breaks = fill_neg(lineas), 
               color = "red4", size = s,  alpha = 0.8) +
  geom_dl(aes(label = ..level.., z = sp), color = "red4", 
          stat = "contour", method = "bottom.pieces", 
          breaks = jumpby(fill_neg(lineas), 2)) +
  
  geom_contour(aes(z = nc), breaks = fill_neg(lineas), 
               color = "blue4", size = s,  alpha = 0.8) +
  geom_dl(aes(label = ..level.., z = nc), color = "blue4", 
          stat = "contour", method = "top.pieces", 
          breaks = jumpby(fill_neg(lineas), 2)) + 
  theme_elio +
  scale_x_continuous(breaks = 1:12, labels = month.abb)

Temperatura

g_data <- rep_lon(sp_nc[variable == "atemp", .(sp = mean(sp), nc = mean(nc)),
                        by = .(lat, lon, lev, estacion)])
maplevs <- plotlevs[plotlevs != 850]
m_data <- g_data[lev %in% maplevs]
M <- m_data[, (max(c(abs(min(nc)), max(nc))))]
M <- ceiling(M)

s <- 0.3
linea <- 2
lineas <- seq(linea, M, by = linea)
# m_data <- m_data[lev == 200]
ggplot(m_data, aes(lon, lat)) +
  # geom_tile(aes(fill = nc)) +
  facet_grid(lev~estacion) +
  HS3+
  geom_contour(aes(z = sp, linetype = as.factor(sign(..level..))),
               breaks = fill_neg(lineas), color = "red4", size = s, alpha = 0.8) +
  # geom_dl(aes(z = sp, label = ..level..), method = "top.pieces", stat = "contour", breaks = fill_neg(jumpby(lineas, 2)), color = "red4") + 
  
  geom_contour(aes(z = nc, linetype = as.factor(sign(..level..))),
               breaks = fill_neg(lineas), color = "blue4", size = s,  alpha = 0.8) +
  # geom_dl(aes(z = nc, label = ..level..), method = "bottom.pieces", stat = "contour", breaks = fill_neg(jumpby(lineas, 2)), color = "blue4") + 
  
  scale_fill_gradient2(high = muted("red"), low = muted("blue"), 
                       name = "NCEP") +
  scale_linetype_manual(values = c(2, 1, 1), guide = "none") +
  coord_map("stereographic", orientation = c(-90,0, 60),
            ylim = c(-90, 0)) +
  scale_x_continuous(breaks = seq(0, 330, by = 30)) +
  # annotate(geom = "text", x = seq(0, 330, by = 30), y = -5, 
  #          label = seq(0, 330, by = 30), color = "gray45") +
  # scale_y_continuous(breaks = seq(0+15, -90+15, by = -15)) +
  # annotate(geom = "text", x = 0, y = seq(0+15, -90+15, by = -15), 
  #          label = seq(0+15, -90+15, by = -15), color = "gray45") +
  theme_bw() + 
  theme(legend.position = "bottom", legend.key.height = unit(5, "points"),
        axis.title = element_blank(), axis.text = element_blank()) +
  ggtitle(paste0("T* para NCEP (azul) y SPEEDY (rojo) \n lineas = ", linea))

linea <- 2
lineas <- seq(linea, M, by = linea)

ggplot(g_data[lat %~% -65], aes(lon, lev)) +
  geom_contour(aes(z = sp, linetype = as.factor(sign(..level..))),
               breaks = fill_neg(lineas), color = "red4", size = s,  alpha = 0.8) +
  geom_contour(aes(z = nc, linetype = as.factor(sign(..level..))), 
               breaks = fill_neg(lineas), color = "blue4", size = s,  alpha = 0.8) +
  scale_linetype_manual(values = c(2, 1, 1), guide = "none") +
  facet_grid(~estacion) +
  scale_y_continuous(trans = "reverselog", breaks = levs) +
  theme_bw() +
  ggtitle(paste0("T* para NCEP (azul) y SPEEDY (rojo) en -65º \n lineas = ", linea))

g_data <- sp_nc_sd[lev %in% plotlevs & variable == "atemp"]

M <- g_data[, (max(c(abs(min(nc)), max(nc))))]
M <- ceiling(M)
s <- 0.7
linea <- 1
lineas <- seq(linea, M, by = linea)

ggplot(g_data, aes(as.numeric(month), lat)) +
  facet_grid(lev~.) +
  geom_contour(aes(z = sp), breaks = fill_neg(lineas), 
               color = "red4", size = s,  alpha = 0.8) +
  geom_dl(aes(label = ..level.., z = sp), color = "red4", 
          stat = "contour", method = "bottom.pieces", 
          breaks = jumpby(fill_neg(lineas), 2)) +
  
  geom_contour(aes(z = nc), breaks = fill_neg(lineas), 
               color = "blue4", size = s,  alpha = 0.8) +
  geom_dl(aes(label = ..level.., z = nc), color = "blue4", 
          stat = "contour", method = "top.pieces", 
          breaks = jumpby(fill_neg(lineas), 2)) + 
  theme_elio + 
  scale_x_continuous(breaks = 1:12, labels = month.abb)

Ondas cuasiestacionarias

colnames(m_ncep.long)[6] <- "value"
colnames(m_speedy.long)[6] <- "value"
m_ncep.long[, modelo := "nc"]
##             lon     lat  lev month variable       value modelo
##       1:   0.00 -87.159 1000   Ene       gh  0.31940276     nc
##       2:   3.75 -87.159 1000   Ene       gh  2.00924682     nc
##       3:   7.50 -87.159 1000   Ene       gh  3.67801808     nc
##       4:  11.25 -87.159 1000   Ene       gh  5.27384569     nc
##       5:  15.00 -87.159 1000   Ene       gh  6.82260795     nc
##      ---                                                      
## 4700156: 341.25  -1.856   10   Dic    atemp  0.19750639     nc
## 4700157: 345.00  -1.856   10   Dic    atemp  0.31581262     nc
## 4700158: 348.75  -1.856   10   Dic    atemp  0.23959038     nc
## 4700159: 352.50  -1.856   10   Dic    atemp  0.22478026     nc
## 4700160: 356.25  -1.856   10   Dic    atemp -0.05261632     nc
m_speedy.long[, modelo := "sp"]
##             lon     lat lev month variable        value modelo
##       1:   0.00 -87.159 925   Ene       gh 736.87003174     sp
##       2:   3.75 -87.159 925   Ene       gh 738.27917887     sp
##       3:   7.50 -87.159 925   Ene       gh 739.81427002     sp
##       4:  11.25 -87.159 925   Ene       gh 741.47501424     sp
##       5:  15.00 -87.159 925   Ene       gh 743.25852458     sp
##      ---                                                      
## 2211836: 341.25  -1.856  30   Dic    atemp  -0.83302260     sp
## 2211837: 345.00  -1.856  30   Dic    atemp  -0.59670922     sp
## 2211838: 348.75  -1.856  30   Dic    atemp  -0.31396132     sp
## 2211839: 352.50  -1.856  30   Dic    atemp   0.04871619     sp
## 2211840: 356.25  -1.856  30   Dic    atemp   0.28689470     sp
sp_nc.long <- rbind(m_ncep.long, m_speedy.long)
remove(m_ncep.long, m_speedy.long)

temp <- sp_nc.long[stringi::stri_sub(variable, 1, 1) == "a"]
qs <- temp[, qs_fit(value, n = 1:4), by = .(lat, lev, month, modelo, variable)]
qs_labels <- c("1" = "QS 1",
               "2" = "QS 2",
               "3" = "QS 3",
               "4" = "QS 4")
mod_labels <- c("nc" = "NCEP", 
                "sp" = "SPEEDY")

R^2

qs_gh <- qs[variable == "agh"]

qs_gh_interp <- qs_gh[!is.na(R.sqr), 
                    interp.dt(lat, lev, R.sqr, linear = T,
                              yo = interp.levs),
                    by = .(month, modelo, variable, QS)]
colnames(qs_gh_interp)[5:7] <- c("lat", "lev", "R.sqr")

temp <- qs_gh[!is.na(Amplitud), 
                    interp.dt(lat, lev, Amplitud, linear = T,
                              yo = interp.levs),
                    by = .(month, modelo, variable, QS)]
colnames(temp)[5:7] <- c("lat", "lev", "Amplitud")

qs_gh_interp$Amplitud <- temp$Amplitud

g_data <- qs_gh_interp
# g_data <- g_data[lat >= -50 & lat < -35 & lev > 100 & QS == 3]
g_nc <- g_data[modelo == "nc"]
g_sp <- g_data[modelo == "sp"]
lineas_width <- 0.25
lineas <- seq(0, 1, by = lineas_width)
s <- 0.7

ggplot(g_nc, aes(lat, lev)) + scale_y_continuous(trans = "reverselog") +
  facet_grid(month ~ QS, labeller = labeller(QS = qs_labels)) +
  geom_tile(aes(fill = cut_width(R.sqr, lineas_width, boundary = 0))) +
  scale_fill_brewer(palette = "Blues", name = "NCEP") +
  geom_contour(aes(z = R.sqr), breaks = lineas, 
               color = "black", data = g_sp, size = s, alpha = 0.7) +
  geom_dl(aes(z = R.sqr, label = ..level..), data = g_sp, color = "black", 
          stat = "contour", method = list("top.pieces", cex = 0.8), 
          breaks = jumpby(lineas, 2)) +
  theme_elio + ggtitle(paste0("R^2 de cada número de onda por latitud y mes \nNCEP (sombreado) y SPEEDY (contornos) - lineas = ", lineas_width))

## Note: no visible binding for global variable 'x' 
## Note: no visible binding for global variable 'width' 
## Note: no visible binding for global variable 'x' 
## Note: no visible binding for global variable 'width' 
## Note: no visible binding for global variable 'y' 
## Note: no visible binding for global variable 'height' 
## Note: no visible binding for global variable 'y' 
## Note: no visible binding for global variable 'height'

Amplitud

s <- 0.7

g_data[lev >= 30 & lev <= 925, MaxAmpl := max(Amplitud), by = .(QS, modelo)]
##         month modelo variable QS        lat  lev     R.sqr Amplitud
##      1:   Ene     nc      agh  1 -87.159000 1000 0.9773782 23.03031
##      2:   Ene     nc      agh  1 -84.971744 1000 0.9283385 32.59569
##      3:   Ene     nc      agh  1 -82.784487 1000 0.8739121 40.36194
##      4:   Ene     nc      agh  1 -80.597231 1000 0.8079078 44.26125
##      5:   Ene     nc      agh  1 -78.409974 1000 0.7822487 45.99572
##     ---                                                            
## 230396:   Dic     sp      agh  4 -10.605026   10        NA       NA
## 230397:   Dic     sp      agh  4  -8.417769   10        NA       NA
## 230398:   Dic     sp      agh  4  -6.230513   10        NA       NA
## 230399:   Dic     sp      agh  4  -4.043256   10        NA       NA
## 230400:   Dic     sp      agh  4  -1.856000   10        NA       NA
##         MaxAmpl
##      1:      NA
##      2:      NA
##      3:      NA
##      4:      NA
##      5:      NA
##     ---        
## 230396:      NA
## 230397:      NA
## 230398:      NA
## 230399:      NA
## 230400:      NA
g_data[, MaxAmpl := max(MaxAmpl, na.rm = T), by = .(QS, modelo)]
##         month modelo variable QS        lat  lev     R.sqr Amplitud
##      1:   Ene     nc      agh  1 -87.159000 1000 0.9773782 23.03031
##      2:   Ene     nc      agh  1 -84.971744 1000 0.9283385 32.59569
##      3:   Ene     nc      agh  1 -82.784487 1000 0.8739121 40.36194
##      4:   Ene     nc      agh  1 -80.597231 1000 0.8079078 44.26125
##      5:   Ene     nc      agh  1 -78.409974 1000 0.7822487 45.99572
##     ---                                                            
## 230396:   Dic     sp      agh  4 -10.605026   10        NA       NA
## 230397:   Dic     sp      agh  4  -8.417769   10        NA       NA
## 230398:   Dic     sp      agh  4  -6.230513   10        NA       NA
## 230399:   Dic     sp      agh  4  -4.043256   10        NA       NA
## 230400:   Dic     sp      agh  4  -1.856000   10        NA       NA
##           MaxAmpl
##      1: 548.81118
##      2: 548.81118
##      3: 548.81118
##      4: 548.81118
##      5: 548.81118
##     ---          
## 230396:  23.22319
## 230397:  23.22319
## 230398:  23.22319
## 230399:  23.22319
## 230400:  23.22319
g_data[, Norm_Ampl := Amplitud/MaxAmpl]
##         month modelo variable QS        lat  lev     R.sqr Amplitud
##      1:   Ene     nc      agh  1 -87.159000 1000 0.9773782 23.03031
##      2:   Ene     nc      agh  1 -84.971744 1000 0.9283385 32.59569
##      3:   Ene     nc      agh  1 -82.784487 1000 0.8739121 40.36194
##      4:   Ene     nc      agh  1 -80.597231 1000 0.8079078 44.26125
##      5:   Ene     nc      agh  1 -78.409974 1000 0.7822487 45.99572
##     ---                                                            
## 230396:   Dic     sp      agh  4 -10.605026   10        NA       NA
## 230397:   Dic     sp      agh  4  -8.417769   10        NA       NA
## 230398:   Dic     sp      agh  4  -6.230513   10        NA       NA
## 230399:   Dic     sp      agh  4  -4.043256   10        NA       NA
## 230400:   Dic     sp      agh  4  -1.856000   10        NA       NA
##           MaxAmpl  Norm_Ampl
##      1: 548.81118 0.04196399
##      2: 548.81118 0.05939327
##      3: 548.81118 0.07354432
##      4: 548.81118 0.08064932
##      5: 548.81118 0.08380974
##     ---                     
## 230396:  23.22319         NA
## 230397:  23.22319         NA
## 230398:  23.22319         NA
## 230399:  23.22319         NA
## 230400:  23.22319         NA
M <- ceiling(g_data[, max(Norm_Ampl, na.rm = T)])
lineas_width <- 0.25
lineas <- c(seq(0, 1, by = lineas_width), 2)

region_3 <- c(latmin = -65, latmax = -40, levmin = 30, levmax = 925)

ggplot(g_data[modelo == "sp"], aes(lat, lev)) + 
                           scale_y_continuous(trans = "reverselog") +
  scale_x_reverse() +
  facet_grid(month ~ QS, labeller = labeller(QS = qs_labels)) +
  
  geom_raster(data = g_data[modelo == "nc"], 
              aes(fill = cut(Norm_Ampl, breaks = lineas))) +

  geom_contour(aes(z = Norm_Ampl), binwidth = lineas_width,
               color = "black", size = s, alpha = 0.7)  + 
  # geom_dl(aes(z = Norm_Ampl, label = ..level..), 
          # binwidth = lineas_width, stat = "contour", method = "top.pieces") +
  annotate(geom = "rect",
           xmin = region_3[1], xmax = region_3[2],
           ymin = region_3[3], ymax = region_3[4],
           fill = NA, color = "black") +
  scale_fill_brewer(palette = "Blues", name = "NCEP") +

  theme_elio + ggtitle(paste0("Amplitud normalizada para cada número de onda por latitud y mes \nNCEP (sombreado) y SPEEDY (contornos) - lineas = ", lineas_width))

QS3

Nos vamos a quedar con la onda 3 y vamos a analizar la zona entre 65º y 40ºS y entre 925 y 30 hPa

g_data <- qs_gh_interp[QS == 3]
region_3 <- c(latmin = -65, latmax = -40, levmin = 30, levmax = 925)
M <- ceiling(g_data[, max(Amplitud, na.rm = T)])
lineas_width <- 10
lineas <- seq(0, M, by = lineas_width)

ggplot(g_data[modelo == "sp"], aes(lat, lev)) + 
  scale_y_continuous(trans = "reverselog") +
  scale_x_reverse() +
  facet_grid(month ~ .) +
  
  geom_raster(data = g_data[modelo == "nc"], 
              aes(fill = cut(Amplitud, breaks = lineas))) +
  
  geom_contour(aes(z = Amplitud), binwidth = lineas_width,
               color = "black", size = s, alpha = 0.7)  + 
  # geom_dl(aes(z = Amplitud, label = ..level..), 
  # binwidth = lineas_width, stat = "contour", method = "top.pieces") +
  annotate(geom = "rect",
           xmin = region_3[1], xmax = region_3[2],
           ymin = region_3[3], ymax = region_3[4],
           fill = NA, color = "black") +
  scale_fill_brewer(palette = "Blues", name = "NCEP") +
  
  theme_elio + ggtitle(paste0("Amplitud de QS3 por latitud y mes \nNCEP (sombreado) y SPEEDY (contornos) - lineas = ", lineas_width))

En esa caja puedo tomar el valor promedio (que sería equivalente a la integral, ya que el área se mantiene igual) o el máximo. El ciclo anual de ambas variables tanto para NCEP como para SPEEDY se muestra en la siguiente figura:

qs3_gh <- qs_gh[QS == 3, ]
qs3_gh[, region := (lat >= region_3[1] & lat <= region_3[2] &
                  lev >= region_3[3] & lev <= region_3[4])]
##           lat  lev month modelo variable  Amplitud       Fase       R.sqr
##    1: -87.159 1000   Ene     nc      agh  1.122137 -1.8936592 0.002320359
##    2: -83.479 1000   Ene     nc      agh  8.521985 -1.6079967 0.042457977
##    3: -79.777 1000   Ene     nc      agh 20.398952 -1.2987663 0.155878052
##    4: -76.070 1000   Ene     nc      agh 23.695890 -1.0897797 0.203066383
##    5: -72.362 1000   Ene     nc      agh 11.180591 -1.0035023 0.080182527
##   ---                                                                    
## 7196: -16.700   30   Dic     sp      agh  4.062847  0.7509012 0.067805701
## 7197: -12.989   30   Dic     sp      agh  5.123547  1.2430367 0.130525363
## 7198:  -9.278   30   Dic     sp      agh  6.235376  1.3842273 0.261522205
## 7199:  -5.567   30   Dic     sp      agh  7.580177  1.3940119 0.363977906
## 7200:  -1.856   30   Dic     sp      agh  9.331789  1.4057353 0.362829617
##       QS region
##    1:  3  FALSE
##    2:  3  FALSE
##    3:  3  FALSE
##    4:  3  FALSE
##    5:  3  FALSE
##   ---          
## 7196:  3  FALSE
## 7197:  3  FALSE
## 7198:  3  FALSE
## 7199:  3  FALSE
## 7200:  3  FALSE
g_data <- qs3_gh[region == T, 
       .(Mean_Ampl = mean(Amplitud), Max_Ampl = max(Amplitud)), 
       by = .(month, modelo)]

g_data[, Area_equiv := Mean_Ampl/Max_Ampl]
##     month modelo Mean_Ampl Max_Ampl Area_equiv
##  1:   Ene     nc 20.255716 40.10602  0.5050542
##  2:   Feb     nc 27.311712 54.67475  0.4995306
##  3:   Mar     nc 22.386954 43.29004  0.5171387
##  4:   Abr     nc 21.981106 40.32700  0.5450717
##  5:   May     nc 17.497137 29.57796  0.5915600
##  6:   Jun     nc 21.332002 36.20745  0.5891606
##  7:   Jul     nc 25.663798 40.75579  0.6296970
##  8:   Ago     nc 27.150506 41.26745  0.6579157
##  9:   Sep     nc 22.962262 36.94286  0.6215616
## 10:   Oct     nc 15.597001 27.31040  0.5711011
## 11:   Nov     nc  6.285683 13.24662  0.4745122
## 12:   Dic     nc  8.606986 18.03878  0.4771379
## 13:   Ene     sp 12.438160 24.64187  0.5047572
## 14:   Feb     sp 13.158030 24.46971  0.5377273
## 15:   Mar     sp 10.394263 18.50236  0.5617804
## 16:   Abr     sp 12.507457 21.28683  0.5875677
## 17:   May     sp 17.952138 30.12335  0.5959543
## 18:   Jun     sp 25.581265 41.36971  0.6183574
## 19:   Jul     sp 17.528378 27.78134  0.6309408
## 20:   Ago     sp 15.358337 20.92206  0.7340740
## 21:   Sep     sp 18.168745 31.32297  0.5800454
## 22:   Oct     sp 17.962713 31.59355  0.5685563
## 23:   Nov     sp 22.239114 38.35022  0.5798953
## 24:   Dic     sp 17.636349 31.80561  0.5545044
##     month modelo Mean_Ampl Max_Ampl Area_equiv
g <- ggplot(g_data, aes(as.numeric(month), Area_equiv)) + 
  geom_line(aes(color = modelo)) +
  scale_color_brewer(type = "qual", palette = 6, labels = c("NCEP", "SPEEDY"),
                     name = "Modelo", direction = -1) +
  scale_x_continuous(breaks = 1:12, labels = month.abb, name = "Mes") +
  scale_y_continuous(limits = c(NA, 1)) +
  ylab("Amplitud media / Amplitud máxima") + 
  theme_elio

g_data <- melt(g_data[, !"Area_equiv", with = F], 
               id.vars = c("modelo", "month"),
               variable.name = "Tipo", value.name = "Amplitud")

ggplot(g_data, aes(as.numeric(month), Amplitud)) +
  geom_line(aes(color = modelo, linetype = Tipo)) +
  scale_x_continuous(breaks = 1:12, labels = month.abb, name = "Mes") +
  scale_color_brewer(type = "qual", palette = 6, labels = mod_labels,
                     name = "Modelo", direction = -1) +
  theme_elio

Se observa que no hay mucha correlación entre los modelos. Mientras que NCEP tiene un ciclo semianual con máximos en febrero y agosto, Speedy tiene los máximos en junio y noviembre. Casi una anticorrelación perfecta.

Otra medida de posible relevancia es la relación entre la amplitud máxima y la media. Ésto da una idea de cuan concentrado está la anomalía de geopotencial. Valores cercanos a 1 implican que la misma está distribuida de forma pareja en todo el área de estudio, mientras que valores menores implican una distribución más concentrada.

g

Esta variable muestra un comportamiento muy parecido en ambos modelos, con un máximo en el invierno y un mínimo en verano, indicando que el patrón de QS3 se encuentra mas localizado en verano que en invierno.

Una cosa que hay que tener cuidado es cómo caracterizar la amplitud de la onda 3 en toda la región mediante un número para hacer una serie temporal. La elección más naive es la media, pero ésto sólo es útil si la variable tiene una distribución simétrica y unimodal. Una variable como la amplitud está acotado inferiormente por el cero, por lo que la primera propiedad no se cumple a priori (aunque puede ser aproximada si la amplitud es grande).

La funciones de distribución para cada mes estimadas a partir de los datos muestran que en varios meses tampoco se cumple la segunda propiedad.

ggplot(qs3_gh[region == T, ], aes(Amplitud)) + 
  geom_density(aes(color = modelo)) + facet_grid(month~.) +
  scale_color_brewer(type = "qual", palette = 6, labels = mod_labels,
                     name = "Modelo", direction = -1) + 
  xlab("Densidad") +
  theme_elio

## Note: no visible binding for global variable 'y' 
## Note: no visible binding for global variable 'x' 
## Note: no visible binding for global variable 'x' 
## Note: no visible binding for global variable 'ymax' 
## Note: no visible binding for global variable 'ymin'

En particular, agosto y noviembre muestran distribuciones multimodales en speedy y ncep respectivamente. Esto quizás podría poner en duda la validez de usar la media como representativo de la actividad de la onda 3 en toda la región.

Serie temporal onda 3

Independientemente de esto, armo una serie temporal mensual con la amplitud de la QS3 para el período 1985-2015 usando datos de NCEP.

# load("NCEP/todo_ncep.rda")
# 
# ncep_r3 <- ncep[lat >= region_3[1] & lat <= region_3[2] &
#                   lev >= region_3[3] & lev <= region_3[4]]
load("NCEP/ncep_r3.rda")
ncep_r3[, agh := anom(gh), by = .(lat, lev, month)]
##             lon     lat lev       time         gh       temp         u
##       1:   0.00 -64.942 925 1985-01-01   531.9597  -4.524887 -1.828449
##       2:   3.75 -64.942 925 1985-01-01   521.6587  -4.454431 -2.684299
##       3:   7.50 -64.942 925 1985-01-01   512.3531  -4.373037 -3.581535
##       4:  11.25 -64.942 925 1985-01-01   505.1609  -4.290367 -4.418997
##       5:  15.00 -64.942 925 1985-01-01   500.0706  -4.237465 -5.100862
##      ---                                                              
## 3499772: 341.25 -42.678  30 2015-12-01 23937.4658 -51.288152  3.380353
## 3499773: 345.00 -42.678  30 2015-12-01 23935.6581 -51.204924  3.560365
## 3499774: 348.75 -42.678  30 2015-12-01 23935.0648 -51.132407  3.708036
## 3499775: 352.50 -42.678  30 2015-12-01 23934.8835 -50.978878  3.867348
## 3499776: 356.25 -42.678  30 2015-12-01 23934.4127 -50.824664  4.091359
##                  v year month        agh
##       1: 4.5704375 1985     1  16.855892
##       2: 4.0168368 1985     1   6.554825
##       3: 3.1929850 1985     1  -2.750730
##       4: 2.2034509 1985     1  -9.942983
##       5: 1.2306730 1985     1 -15.033255
##      ---                                
## 3499772: 0.8035898 2015    12  -4.778157
## 3499773: 0.6716072 2015    12  -6.585878
## 3499774: 0.4757195 2015    12  -7.179177
## 3499775: 0.2354143 2015    12  -7.360489
## 3499776: 0.1725060 2015    12  -7.831237
qs3 <- ncep_r3[, qs_fit(agh, n = 3), by = .(lat, lev, time)]
qs3[, month := as.numeric(stringi::stri_sub(time, 6, 7))]
##            lat lev       time  Amplitud         Fase       R.sqr QS month
##     1: -64.942 925 1985-01-01 11.803611  2.074606069 0.115521622  3     1
##     2: -61.232 925 1985-01-01 17.083049  2.266734636 0.114262180  3     1
##     3: -57.521 925 1985-01-01 22.758725  2.442488409 0.142396274  3     1
##     4: -53.810 925 1985-01-01 24.283161  2.554906552 0.155590425  3     1
##     5: -50.099 925 1985-01-01 21.845354  2.589264691 0.154252904  3     1
##    ---                                                                   
## 36452: -57.521  30 2015-12-01 25.537137 -0.337952733 0.012160818  3    12
## 36453: -53.810  30 2015-12-01 20.496370 -0.387828235 0.011492796  3    12
## 36454: -50.099  30 2015-12-01 15.693286 -0.406047472 0.011419940  3    12
## 36455: -46.389  30 2015-12-01 10.640019 -0.271296775 0.010503649  3    12
## 36456: -42.678  30 2015-12-01  5.997953  0.006903003 0.008447028  3    12
qs3[, month := factor(month, levels = c(12, 1:11), ordered = T)]
##            lat lev       time  Amplitud         Fase       R.sqr QS month
##     1: -64.942 925 1985-01-01 11.803611  2.074606069 0.115521622  3     1
##     2: -61.232 925 1985-01-01 17.083049  2.266734636 0.114262180  3     1
##     3: -57.521 925 1985-01-01 22.758725  2.442488409 0.142396274  3     1
##     4: -53.810 925 1985-01-01 24.283161  2.554906552 0.155590425  3     1
##     5: -50.099 925 1985-01-01 21.845354  2.589264691 0.154252904  3     1
##    ---                                                                   
## 36452: -57.521  30 2015-12-01 25.537137 -0.337952733 0.012160818  3    12
## 36453: -53.810  30 2015-12-01 20.496370 -0.387828235 0.011492796  3    12
## 36454: -50.099  30 2015-12-01 15.693286 -0.406047472 0.011419940  3    12
## 36455: -46.389  30 2015-12-01 10.640019 -0.271296775 0.010503649  3    12
## 36456: -42.678  30 2015-12-01  5.997953  0.006903003 0.008447028  3    12
qs3[, year := as.numeric(stringi::stri_sub(time, 1, 4))]
##            lat lev       time  Amplitud         Fase       R.sqr QS month
##     1: -64.942 925 1985-01-01 11.803611  2.074606069 0.115521622  3     1
##     2: -61.232 925 1985-01-01 17.083049  2.266734636 0.114262180  3     1
##     3: -57.521 925 1985-01-01 22.758725  2.442488409 0.142396274  3     1
##     4: -53.810 925 1985-01-01 24.283161  2.554906552 0.155590425  3     1
##     5: -50.099 925 1985-01-01 21.845354  2.589264691 0.154252904  3     1
##    ---                                                                   
## 36452: -57.521  30 2015-12-01 25.537137 -0.337952733 0.012160818  3    12
## 36453: -53.810  30 2015-12-01 20.496370 -0.387828235 0.011492796  3    12
## 36454: -50.099  30 2015-12-01 15.693286 -0.406047472 0.011419940  3    12
## 36455: -46.389  30 2015-12-01 10.640019 -0.271296775 0.010503649  3    12
## 36456: -42.678  30 2015-12-01  5.997953  0.006903003 0.008447028  3    12
##        year
##     1: 1985
##     2: 1985
##     3: 1985
##     4: 1985
##     5: 1985
##    ---     
## 36452: 2015
## 36453: 2015
## 36454: 2015
## 36455: 2015
## 36456: 2015

Meses de verano donde activo e inactivo se definen si la amplitud es mayor, menor o que la media (del mes) en 1 desvío estandar (del mes):

qs3_mean <- qs3[, .(Mean_Ampl = mean(Amplitud),
                    Max_Ampl = max(Amplitud)), 
                by = .(month, year, time)]

qs3_mean[, `:=`(MeanMean_Ampl = mean(Mean_Ampl),
                SDMean_Ampl = sd(Mean_Ampl),
                MeanMax_Ampl = mean(Max_Ampl),
                SDMax_Ampl = sd(Max_Ampl)),
                by = .(month)]
##      month year       time Mean_Ampl  Max_Ampl MeanMean_Ampl SDMean_Ampl
##   1:     1 1985 1985-01-01  24.52334  58.86156      28.67268    11.79553
##   2:     2 1985 1985-02-01  30.08895  67.91443      34.54357    12.04654
##   3:     3 1985 1985-03-01  36.82800  64.30854      33.86206    12.32747
##   4:     4 1985 1985-04-01  38.25564  72.74778      37.35963    17.38566
##   5:     5 1985 1985-05-01  21.27393  32.00394      36.03709    16.72060
##  ---                                                                    
## 368:     8 2015 2015-08-01  24.59589  39.53484      50.04477    23.54863
## 369:     9 2015 2015-09-01  69.22008 100.90178      44.63204    24.43700
## 370:    10 2015 2015-10-01  44.67534  64.20312      34.86898    16.19144
## 371:    11 2015 2015-11-01  15.24799  35.95592      34.99821    15.94231
## 372:    12 2015 2015-12-01  15.31068  28.22762      26.72547    15.62461
##      MeanMax_Ampl SDMax_Ampl
##   1:     59.79659   22.14225
##   2:     69.93723   23.27483
##   3:     64.56587   23.92974
##   4:     68.57183   28.63779
##   5:     61.30037   23.99550
##  ---                        
## 368:     75.21750   33.08701
## 369:     71.29075   34.84362
## 370:     59.18704   20.95578
## 371:     60.07507   23.82482
## 372:     51.54154   26.40450
shift_season <- function(month, year) {
  ifelse(month == 12, year + 1, year)
}

qs3_mean[, Año_estacion := shift_season(month, year)]
##      month year       time Mean_Ampl  Max_Ampl MeanMean_Ampl SDMean_Ampl
##   1:     1 1985 1985-01-01  24.52334  58.86156      28.67268    11.79553
##   2:     2 1985 1985-02-01  30.08895  67.91443      34.54357    12.04654
##   3:     3 1985 1985-03-01  36.82800  64.30854      33.86206    12.32747
##   4:     4 1985 1985-04-01  38.25564  72.74778      37.35963    17.38566
##   5:     5 1985 1985-05-01  21.27393  32.00394      36.03709    16.72060
##  ---                                                                    
## 368:     8 2015 2015-08-01  24.59589  39.53484      50.04477    23.54863
## 369:     9 2015 2015-09-01  69.22008 100.90178      44.63204    24.43700
## 370:    10 2015 2015-10-01  44.67534  64.20312      34.86898    16.19144
## 371:    11 2015 2015-11-01  15.24799  35.95592      34.99821    15.94231
## 372:    12 2015 2015-12-01  15.31068  28.22762      26.72547    15.62461
##      MeanMax_Ampl SDMax_Ampl Año_estacion
##   1:     59.79659   22.14225         1985
##   2:     69.93723   23.27483         1985
##   3:     64.56587   23.92974         1985
##   4:     68.57183   28.63779         1985
##   5:     61.30037   23.99550         1985
##  ---                                     
## 368:     75.21750   33.08701         2015
## 369:     71.29075   34.84362         2015
## 370:     59.18704   20.95578         2015
## 371:     60.07507   23.82482         2015
## 372:     51.54154   26.40450         2016
qs3_mean[, Anom := ifelse(Mean_Ampl > MeanMean_Ampl + SDMean_Ampl,
                          "Activo", 
                          ifelse(Mean_Ampl < MeanMean_Ampl - SDMean_Ampl, 
                                 "Inactivo", "Normal"))]
##      month year       time Mean_Ampl  Max_Ampl MeanMean_Ampl SDMean_Ampl
##   1:     1 1985 1985-01-01  24.52334  58.86156      28.67268    11.79553
##   2:     2 1985 1985-02-01  30.08895  67.91443      34.54357    12.04654
##   3:     3 1985 1985-03-01  36.82800  64.30854      33.86206    12.32747
##   4:     4 1985 1985-04-01  38.25564  72.74778      37.35963    17.38566
##   5:     5 1985 1985-05-01  21.27393  32.00394      36.03709    16.72060
##  ---                                                                    
## 368:     8 2015 2015-08-01  24.59589  39.53484      50.04477    23.54863
## 369:     9 2015 2015-09-01  69.22008 100.90178      44.63204    24.43700
## 370:    10 2015 2015-10-01  44.67534  64.20312      34.86898    16.19144
## 371:    11 2015 2015-11-01  15.24799  35.95592      34.99821    15.94231
## 372:    12 2015 2015-12-01  15.31068  28.22762      26.72547    15.62461
##      MeanMax_Ampl SDMax_Ampl Año_estacion     Anom
##   1:     59.79659   22.14225         1985   Normal
##   2:     69.93723   23.27483         1985   Normal
##   3:     64.56587   23.92974         1985   Normal
##   4:     68.57183   28.63779         1985   Normal
##   5:     61.30037   23.99550         1985   Normal
##  ---                                              
## 368:     75.21750   33.08701         2015 Inactivo
## 369:     71.29075   34.84362         2015   Activo
## 370:     59.18704   20.95578         2015   Normal
## 371:     60.07507   23.82482         2015 Inactivo
## 372:     51.54154   26.40450         2016   Normal
# Alternativa: usando terciles. 
# qs3_mean[, Anom := cut(Mean_Ampl, 
#                breaks = quantile(Mean_Ampl, seq(0, 1, length.out = 4)),
#                include.lowest = T, labels = c("Inactivo", "Normal", "Activo")),
#          by = month]

years <- seq(1984, 2016, by = 1)
lab <- stri_sub(years, 3, 4)
labs <- paste0(lab, "/", shift(lab, type = "lead"))
n <- length(years) 
years <- years[-c(1)]
labs <- labs[-c(n)]
# years <- unique(qs3_mean$Año_estacion)


g <- ggplot(qs3_mean, 
       aes(Año_estacion, Mean_Ampl)) + 
  geom_col(aes(fill = Anom)) +
  geom_line(aes(y = MeanMean_Ampl), linetype = 2, color = "gray4") +
  scale_x_continuous(breaks = jumpby(years, 2),
                     labels = jumpby(labs, 2),
                     name = "Año de fin de estación") +
  scale_fill_brewer(palette = "Set1")  + 
  ylab("Amplitud media QS3") +
  theme_elio + 
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))
g + facet_grid(month~., labeller = labeller(month = month.abb))

## Note: no visible binding for global variable 'y' 
## Note: no visible binding for global variable 'y' 
## Note: no visible binding for global variable 'x' 
## Note: no visible binding for global variable 'width' 
## Note: no visible binding for global variable 'x' 
## Note: no visible binding for global variable 'width'

La gran mayoría de los meses de veranos son normales con 5-6 veranos activo o inactivos.

Se observa que no hay muchos meses con amplitud mucho mayor que la media y que no están agrupados en un mismo verano. Es decir, no hay mucha correlación entre la amplitud de un mes y el siguiente.

Para ver esto último con más exactitud, se puede estimar la función de autocorrelación.

autocor <- acf.sig(qs3_mean$Mean_Ampl, method = "large.lag")

ggplot(autocor, aes(lag, acf)) + 
  geom_line() + geom_point() +
  geom_line(aes(y = sig.cut), linetype = 2) +
  geom_line(aes(y = -sig.cut), linetype = 2) +
  geom_hline(yintercept = 0) +
  scale_x_continuous(limits = c(0, 12*3),
                     name = "Lag (meses)", 
                     breaks = seq(0, 12*3, by = 6)) +
  ylab("Autocorrelación mensual QS3") +
  theme_elio

Se ve que la autocorrelación a lag 1 no es significativa (utilizando método large lag), pero además se ve un evidente ciclo anual.

ggplot(qs3_mean, aes(factor(month, labels = month.abb), Mean_Ampl)) + 
  geom_boxplot() + 
  xlab("Mes") + ylab("Amplitud media QS3") +
  theme_elio

## Note: no visible binding for global variable 'weight' 
## Note: no visible binding for global variable 'xmin' 
## Note: no visible binding for global variable 'xmax' 
## Note: no visible binding for global variable 'y' 
## Note: no visible binding for global variable 'size'

¿Tiene otros ciclos? Le quito el ciclo anual restando la media anual y hago fourier.

fourier <- as.dt(convert.fft(fft(qs3_mean[, .(Mean_Ampl - MeanMean_Ampl)]$V1)))

ggplot(fourier, aes(freq, ampl)) + 
  geom_line() +
  theme_elio + xlab("Frecuencia") + ylab("Amplitud")

No, no hay ningún ciclo obvio.

qs3 <- merge(qs3, qs3_mean[, .(month, year, Anom)])
# 
# g_data <- qs3[, .(Amplitud = mean(Amplitud)), by = .(lat, lev, Anom, month)]
# 
# g_data <- g_data[, interp.dt(lat, lev, Amplitud, linear = T,
#                               yo = interp.levs),
#                     by = .(Anom, month)]
# colnames(g_data) <- c("Anom", "month", "lat", "lev", "Amplitud")
# # g_data <- g_data[month %in% 1:2]
# ggplot(g_data, aes(lat, lev)) + 
#   scale_y_continuous(trans = "reverselog") +
#   geom_contour(aes(z = Amplitud), binwidth = 10, 
#                color = "black") +
#   geom_dl(aes(z = Amplitud, label = ..level..),
#           stat = "contour", binwidth = 20, color = "black",
#           method = "top.pieces") +
#   facet_grid(month~Anom, labeller = labeller(month = month.abb)) +
#   theme_elio
#   

Esto es la distribución espacial de la amplitud en la región de estudio para cada mes y clasificación. No sé si sirve para mucho. Se ve que en febrero en algunos meses hay más separación que en otros.

ggplot(qs3, aes(Amplitud)) + 
  geom_density(aes(color = Anom)) +
  facet_wrap(~month, ncol = 4, labeller = labeller(month = month.abb)) +
  theme_elio

Composiciones de variables.

Ahora la idea es ver los campos de variables que podrían ser forantes y comparar años activos vs. inactivos.

SST

load("NCEP/sst_todo.rda")

## Agrego máscara de océanos
library(maptools)
library(maps)

make_mask <- function(lat, lon) {
  seamask <- map("world2", fill=TRUE, col = "transparent", plot = F)
  IDs <- sapply(strsplit(seamask$names, ":"), function(x) x[1])
  seamask <- map2SpatialPolygons(seamask, IDs = IDs,
                                 proj4string = CRS("+proj=longlat +datum=WGS84"))
  
  points <- SpatialPoints(expand.grid(lon, lat), 
                          proj4string = CRS(proj4string(seamask)))
  sea <-  is.na(over(points, seamask))
  return(sea)
}

lat <- unique(sst$lat)
lon <- unique(sst$lon)
sea <- make_mask(lat, lon)
invisible(sst[, sea := sea])

sst <- merge(sst, qs3_mean[, .(time, Anom, Mean_Ampl)], all = T)
invisible({
  sst[, month := as.numeric(stringi::stri_sub(time, 6, 7))]
  sst[, month := factor(month, levels = c(12, 1:11), ordered = T)]
  sst[sea == F, sst := NA]
  sst[, msst := mean(sst, na.rm = T), by = .(lat, month)]
  sst[, asst := sst - msst]
})
sst_comp <- sst[!is.na(Anom), 
           .(asst = mean(asst, na.rm = T)), 
           by = .(lat, lon, Anom, month)]
tmp <- sst_comp[, .(Anom = "Diferencia",
                    asst = asst[Anom == "Inactivo"] - asst[Anom == "Activo"]),
                by = .(lat, lon, month)]
sst_comp <- rbind(sst_comp, tmp)
remove(tmp)

gdata <- sst_comp[Anom %in% c("Activo", "Inactivo")]


g <- ggplot(gdata, aes(lon, lat)) + 
  geom_contour(aes(z = asst, color = ..level..), binwidth = 2.5) +
  world3 +
  coord_quickmap() +
  scale_color_gradient2(high = muted("red"), low = muted("blue")) +
  theme_elio

DEF: (líneas = 2.5)

g + facet_grid_paginate(month~Anom, labeller = labeller(month = month.abb), 
                         nrow = 3, ncol = 2,
                         page = 1)

MAM

g + facet_grid_paginate(month~Anom, labeller = labeller(month = month.abb), 
                         nrow = 3, ncol = 2,
                         page = 2)

JJA

g + facet_grid_paginate(month~Anom, labeller = labeller(month = month.abb), 
                         nrow = 3, ncol = 2,
                         page = 3)

SON:

g + facet_grid_paginate(month~Anom, labeller = labeller(month = month.abb), 
                         nrow = 3, ncol = 2,
                         page = 4)

No se puede ver a ojo mucha diferencia. Restando ambos campos (íneas = 0.5):

lineas <- seq(-3, 3, by = 0.5)
gdata <- sst_comp[Anom == "Diferencia"]
ggplot(gdata, aes(lon, lat)) + 
  geom_contour(aes(z = asst, color = ..level..), breaks = lineas) + 
  world3 +
  coord_quickmap() +
  facet_wrap(~month, labeller = labeller(month = month.abb), ncol = 3) + 
  scale_color_gradient2(high = muted("red"), low = muted("blue")) +
  theme_elio

En algunos meses se ve que la zona del pacífico ecuatorial tiene grandes diferencias, con nomviembre, diciemre y enero caracterizado por diferencias positivas (años activos con mayor temperatura que inactivos) aunque lo contrario en septiembre y febrero.

Al hacer la clasificación binaria entre activos e inactivos se pierde información. Una alternativa es hacer la regresión de la SST en la amplitud de la onda QS3. Es decir, modelar \(QS3 = m\times SST + b + \epsilon\) para cada punto del globo.

sst[!is.na(Mean_Ampl) & !is.na(asst), sd := sd(asst), by = .(lat, lon, month)]
##                time    lon     lat       sst   sea   Anom Mean_Ampl month
##       1: 1984-01-01   0.00 -87.159        NA FALSE     NA        NA     1
##       2: 1984-01-01   3.75 -87.159        NA FALSE     NA        NA     1
##       3: 1984-01-01   7.50 -87.159        NA FALSE     NA        NA     1
##       4: 1984-01-01  11.25 -87.159        NA FALSE     NA        NA     1
##       5: 1984-01-01  15.00 -87.159        NA FALSE     NA        NA     1
##      ---                                                                 
## 1769468: 2015-12-01 341.25  87.159 -1.742222  TRUE Normal  15.31068    12
## 1769469: 2015-12-01 345.00  87.159 -1.776820  TRUE Normal  15.31068    12
## 1769470: 2015-12-01 348.75  87.159 -1.790000  TRUE Normal  15.31068    12
## 1769471: 2015-12-01 352.50  87.159 -1.790000  TRUE Normal  15.31068    12
## 1769472: 2015-12-01 356.25  87.159 -1.790000  TRUE Normal  15.31068    12
##               msst         asst           sd
##       1:       NaN           NA           NA
##       2:       NaN           NA           NA
##       3:       NaN           NA           NA
##       4:       NaN           NA           NA
##       5:       NaN           NA           NA
##      ---                                    
## 1769468: -1.786385  0.044162663 0.0140890950
## 1769469: -1.786385  0.009565163 0.0042702613
## 1769470: -1.786385 -0.003614836 0.0000000000
## 1769471: -1.786385 -0.003614836 0.0000000000
## 1769472: -1.786385 -0.003614836 0.0008876992
# sst_lm <-
#   sst[!is.na(Mean_Ampl) & !is.na(asst) & sd != 0,
#       {
#         a <- summary(lm(Mean_Ampl ~ asst))
#         list(estimate  = a$coefficients[2, 1],
#              se        = a$coefficients[2, 2])
#       },
#       by = .(lon, lat, month)]
# save(sst_lm, file = "sst_lm.rda")
load("sst_lm.rda")

# sst_ARMA <-
#   sst[!is.na(Mean_Ampl) & !is.na(asst) & sd != 0,
#       {
#         a <- arima(Mean_Ampl, xreg = asst,
#                    order = c(1, 0, 1), method = "ML")
#         list(estimate = coef(a)[[4]],
#              se       = sqrt(diag(a$var.coef))[4])
#       },
#       by = .(lon, lat, month)]
# save(sst_ARMA, file = "sst_ARMA.rda")
load("sst_ARMA.rda")

sst_ARMA[, model := "ARMA"]
##           lon     lat month      estimate          se model
##     1: 165.00 -76.070     1    -0.5710952    6.160772  ARMA
##     2: 168.75 -76.070     1    -1.0979065    6.955729  ARMA
##     3: 172.50 -76.070     1    16.3959292    6.298035  ARMA
##     4: 176.25 -76.070     1     5.9687953    9.500008  ARMA
##     5: 180.00 -76.070     1    -1.4928169    7.130592  ARMA
##    ---                                                     
## 35617: 333.75  87.159    12  -227.2418053  576.201620  ARMA
## 35618: 337.50  87.159    12  -209.5381702  184.037547  ARMA
## 35619: 341.25  87.159    12  -265.9381561  249.217011  ARMA
## 35620: 345.00  87.159    12  -464.5931503  675.701240  ARMA
## 35621: 356.25  87.159    12 -1662.0564847 2741.928059  ARMA
sst_lm[, model := "lm"]
##           lon     lat month     estimate          se model
##     1: 165.00 -76.070     1     3.465502    7.369890    lm
##     2: 168.75 -76.070     1     6.279632    7.834251    lm
##     3: 172.50 -76.070     1     9.603780    9.051243    lm
##     4: 176.25 -76.070     1    13.800397    9.930859    lm
##     5: 180.00 -76.070     1     7.110357    7.979328    lm
##    ---                                                    
## 35617: 333.75  87.159    12  -418.695999  589.880315    lm
## 35618: 337.50  87.159    12   -66.482428  156.488025    lm
## 35619: 341.25  87.159    12   -67.993298  205.546154    lm
## 35620: 345.00  87.159    12  -288.471974  677.332404    lm
## 35621: 356.25  87.159    12 -1405.224116 3258.036442    lm
sst_regr <- rbind(sst_ARMA, sst_lm)
remove(sst_ARMA, sst_lm)

sst_regr[, t := abs(estimate/se)]
##           lon     lat month      estimate          se model          t
##     1: 165.00 -76.070     1    -0.5710952    6.160772  ARMA 0.09269863
##     2: 168.75 -76.070     1    -1.0979065    6.955729  ARMA 0.15784205
##     3: 172.50 -76.070     1    16.3959292    6.298035  ARMA 2.60334024
##     4: 176.25 -76.070     1     5.9687953    9.500008  ARMA 0.62829372
##     5: 180.00 -76.070     1    -1.4928169    7.130592  ARMA 0.20935386
##    ---                                                                
## 71238: 333.75  87.159    12  -418.6959993  589.880315    lm 0.70979822
## 71239: 337.50  87.159    12   -66.4824280  156.488025    lm 0.42484035
## 71240: 341.25  87.159    12   -67.9932976  205.546154    lm 0.33079333
## 71241: 345.00  87.159    12  -288.4719743  677.332404    lm 0.42589425
## 71242: 356.25  87.159    12 -1405.2241162 3258.036442    lm 0.43131013
ggplot(sst_regr[model == "lm" & abs(estimate) < 40], aes(lon, lat)) + 
  geom_raster(aes(fill = estimate)) +
  # geom_contour(aes(z = abs(estimate/se)), breaks = 2:100, color = "black") +
  geom_point(size = 0.1, shape = 3, alpha = 0.4,
             data = sst_regr[model == "lm" & t > 2]) +
  world3 + coord_quickmap() +
  theme_elio +
  facet_wrap(~month, labeller = labeller(month = month.abb), ncol = 3) + 
  scale_fill_gradient2(high = muted("red"), low = muted("blue"),
                       name = "Regresión (LM)", limits = c(-40, 40))

Si quiero ser un poco más riguroso y tener en cuenta la estructura de correlación que seguro existe en las series, puedo usar un modelo ARMA(1, 1) para estimar la pendiente.

ggplot(sst_regr[model == "ARMA" & abs(estimate) < 40], aes(lon, lat)) + 
  geom_raster(aes(fill = estimate)) +
  # geom_contour(aes(z = abs(estimate/se)), breaks = 2:100, color = "black") +
  geom_point(size = 0.1, shape = 3, alpha = 0.4,
             data = sst_regr[model == "ARMA" & t > 2]) +
  world3 + coord_quickmap() +
  theme_elio +
  facet_wrap(~month, labeller = labeller(month = month.abb), ncol = 3) + 
  scale_fill_gradient2(high = muted("red"), low = muted("blue"),
                       name = "Regresión (ARMA)", limits = c(-40, 40))

En ambas figuras las cruces marcan regiones donde la magnitud de la regresión es mayor a 2 veces el error estandar de la misma (~ p-valor < 0.05).

(Alternativa: hacer test t-student para cada punto de latitud entre años activos e inactivos. )

# test <- 
#   sst[!is.na(Mean_Ampl) & !is.na(asst) & sd != 0,
#       t.test(asst[Anom == "Activo"], y = asst[Anom == "Inactivo"]), 
#       by = .(lat, lon, month)]
# 
# ggplot(test, aes(lon, lat)) +
#   geom_raster(aes(fill = estimate))  +
#   geom_point(aes(alpha = p.value < 0.05), size = 0.1, shape = 3) +
#   scale_alpha_discrete(range = c(0, 0.7)) +
#   scale_fill_gradient2(high = muted("red"), low = muted("blue")) +
#   facet_wrap(~month, labeller = labeller(month = month.abb), ncol = 3) +
#   coord_quickmap() +
#   world3 +
#   theme_elio